geom_smooth(aes(x=GD, y=Alt), method = "lm", color="darkred")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$GD),max(city_data$GD))+
geom_point(aes(x=APG, y=Alt), size=3, color="yellow")+
geom_smooth(aes(x=APG, y=Alt), method = "lm", color="darkgreen")
ggplot(data=city_data)+
geom_point(aes(x=GD, y=Alt), size=3, color="darkblue")+
geom_smooth(aes(x=GD, y=Alt), method = "lm", color="darkred")+
ylim(0,max(city_data$Alt))+
xlim(0,max(city_data$GD))+
geom_point(aes(x=APG, y=Alt), size=3, color="yellow")+
geom_smooth(aes(x=APG, y=Alt), method = "lm", color="darkgreen")
ggplot(data=city_data)+
geom_point(aes(x=GD, y=Alt), size=3, color="darkblue")+
geom_smooth(aes(x=GD, y=Alt), method = "lm", color="darkred")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$GD),max(city_data$GD))+
theme_bw()
ggplot(data=city_data)+
geom_point(aes(x=GD, y=Alt), size=3, color="darkblue")+
geom_smooth(aes(x=GD, y=Alt), method = "lm", color="darkred")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$GD),max(city_data$GD))+
theme_bw() +
ylab("Altitude") + xlab("GD %")+
labs(title = "Districts", caption = "2024 CEC Official results")
ggplot(data=city_data)+
geom_point(aes(x=GD, y=Alt), size=3, color="darkblue")+
geom_smooth(aes(x=GD, y=Alt), method = "lm", color="darkred")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$GD),max(city_data$GD))+
theme_bw() +
ylab("Altitude") + xlab("GD %")+
labs(title = "Districts", subtitle = "2024 CEC Official results")
corrplot::corrplot(city_data$GD)
corrplot::corrplot(city_data)
View(city_data)
corrplot::corrplot(city_data[,c(2,4)])
city_data[,c(2,4)]
cor(city_data$GD, city_data$Alt)
ggplot(data=city_data)+
geom_point(aes(x=GD, y=Alt), size=3, color="#017e84")+
geom_smooth(aes(x=GD, y=Alt), method = "lm", color="#e02626")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$GD),max(city_data$GD))+
theme_bw() +
ylab("Altitude") + xlab("GD %")+
labs(title = "Districts", subtitle = "2024 CEC Official results")
ggplot(data=city_data)+
geom_point(aes(x=GD, y=Alt), size=3, color="#0044d5")+
geom_smooth(aes(x=GD, y=Alt), method = "lm", color="#e02626")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$GD),max(city_data$GD))+
theme_bw() +
ylab("Altitude") + xlab("GD %")+
labs(title = "Districts", subtitle = "2024 CEC Official results")
ggplot(data=city_data)+
geom_point(aes(x=APG, y=Alt), size=3, color="#0044d5")+
geom_smooth(aes(x=APG, y=Alt), method = "lm", color="#e02626")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$APG),max(city_data$APG))+
theme_bw() +
ylab("Altitude") + xlab("APG %")+
labs(title = "Districts", subtitle = "2024 CEC Official results")
cor(city_data$APG, city_data$Alt)
ggplot(data=city_data)+
geom_point(aes(x=APG, y=Alt), size=3, color="#d5b800")+
geom_smooth(aes(x=APG, y=Alt), method = "lm", color="#e02626")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$APG),max(city_data$APG))+
theme_bw() +
ylab("Altitude") + xlab("APG %")+
labs(title = "Districts", subtitle = "2024 CEC Official results")
ggplot(data=city_data)+
geom_point(aes(x=GD, y=Alt), size=3, color="#0044d5")+
geom_smooth(aes(x=GD, y=Alt), method = "lm", color="#e02626")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$GD),max(city_data$GD))+
theme_bw() +
ylab("Altitude") + xlab("GD %")+
labs(title = "Districts", subtitle = "2024 CEC Official results")
ggplot(data=city_data)+
geom_point(aes(x=APG, y=Alt), size=3, color="#d5b800")+
geom_smooth(aes(x=APG, y=Alt), method = "lm", color="#e02626")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$APG),max(city_data$APG))+
theme_bw() +
ylab("Altitude") + xlab("APG %")+
labs(title = "Districts", subtitle = "2024 CEC Official results")
ggplot(data=city_data)+
geom_point(aes(x=GD, y=Alt), size=3, color="#0044d5")+
geom_smooth(aes(x=GD, y=Alt), method = "lm", color="#e02626")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$GD),max(city_data$GD))+
theme_bw() +
ylab("Altitude") + xlab("GD %")+
labs(title = "GD vote share x Altitude", subtitle = "Districts \n2024 CEC Official results")
ggplot(data=city_data)+
geom_point(aes(x=APG, y=Alt), size=3, color="#d5b800")+
geom_smooth(aes(x=APG, y=Alt), method = "lm", color="#e02626")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$APG),max(city_data$APG))+
theme_bw() +
ylab("Altitude") + xlab("APG %")+
labs(title = "APG vote share x Altitude", subtitle = "Districts \n2024 CEC Official results")
ggplot(data=city_data)+
geom_point(aes(x=GD, y=Alt), size=3, color="#0044d5")+
geom_smooth(aes(x=GD, y=Alt), method = "lm", color="#e02626")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$GD),max(city_data$GD))+
theme_bw() +
ylab("Altitude") + xlab("GD %")+
labs(title = "GD vote share x Altitude", subtitle = "Districts \n2024 CEC Official results",
caption = "r=0.63")
ggplot(data=city_data)+
geom_point(aes(x=APG, y=Alt), size=3, color="#d5b800")+
geom_smooth(aes(x=APG, y=Alt), method = "lm", color="#e02626")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$APG),max(city_data$APG))+
theme_bw() +
ylab("Altitude") + xlab("APG %")+
labs(title = "APG vote share x Altitude", subtitle = "Districts \n2024 CEC Official results",
caption = "r=-0.12")
ggplot(data=city_data)+
geom_point(aes(x=GD, y=Alt), size=3, color="#0044d5")+
geom_smooth(aes(x=GD, y=Alt), method = "lm", color="#e02626")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$GD),max(city_data$GD))+
theme_bw() +
ylab("Altitude") + xlab("GD %")+
labs(title = "GD vote share x Altitude", subtitle = "Districts \n2024 CEC Official results",
caption = "r = 0.63")
ggplot(data=city_data)+
geom_point(aes(x=APG, y=Alt), size=3, color="#d5b800")+
geom_smooth(aes(x=APG, y=Alt), method = "lm", color="#e02626")+
ylim(0,max(city_data$Alt))+
xlim(min(city_data$APG),max(city_data$APG))+
theme_bw() +
ylab("Altitude") + xlab("APG %")+
labs(title = "APG vote share x Altitude", subtitle = "Districts \n2024 CEC Official results",
caption = "r = -0.12")
sqrt(100)
sqrt(1034)
sqrt(1024)
#######################################################################
#Name of code file: mobilization_analysis.R
#Data In: st_panel.RDS
#Data Out: Figures 5, A10; Tables 2, A9, A10, A11, A12, A13, A14
######################################################################
# ----------------------------------------------------------------------
# load all relevant packages
# ----------------------------------------------------------------------
## If you don't have packages installed, uncomment and install
## install.packages("tidyverse")
## install.packages("ggplot2")
## install.packages("stargazer")
## install.packages("xtable")
## install.packages("texreg")
## install.packages("readstata13")
## install.packages("estimatr")
## install.packages("naniar")
## install.packages("lfe")
## install.packages("multiwayvcov")
## install.packages("lmtest")
## install.packages("xlsx")
## install.packages("readxl")
## install.packages("lubridate")
## install.packages("panelView")
## install.packages("zoo")
## Note that the Zelig package, which we use for making first differences plots,
## has been removed from CRAN - in order to install, you should first download
## the latest version from the Archive (5.1.7 as of this script) and manually
## install it using the following code - the first line should be the filepath
## to the .tar.gz file that you downloaded, the second line uses this filepath
## to install the package
## zeligpackage <- "~/Downloads/Zelig_5.1.7.tar.gz"
## install.packages(zeligpackage, repos = NULL, type = "source")
## install.packages("MASS")
## install.packages("sandwich")
## install.packages("lmtest")
## install.packages("bucky")
## install.packages("scales")
library("tidyverse")
library("ggplot2")
library("stargazer")
library("xtable")
library("texreg")
library("readstata13")
library("estimatr")
library("naniar")
library("lfe")
library("multiwayvcov")
library("lmtest")
library("xlsx")
library("readxl")
library("lubridate")
library("panelView")
library("zoo")
library("Zelig")
library("MASS")
library("sandwich")
library("lmtest")
# library("bucky")
library("scales")
## On the computer used to generate the results of this script, below is what is
## printed to the console when sessionInfo() is run
## sessionInfo()
## > sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Monterey 12.1
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRlapack.dylib
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base
## other attached packages:
##  [1] scales_1.1.1       bucky_1.0.6        sandwich_3.0-1     MASS_7.3-54
##  [5] Zelig_5.1.7        survival_3.2-13    panelView_1.1.5    lubridate_1.8.0
##  [9] readxl_1.3.1       xlsx_0.6.5         lmtest_0.9-39      zoo_1.8-9
## [13] multiwayvcov_1.2.3 lfe_2.8-7.1        Matrix_1.3-4       naniar_0.6.1
## [17] estimatr_0.30.4    readstata13_0.10.0 texreg_1.37.5      xtable_1.8-4
## [21] stargazer_5.2.2    forcats_0.5.1      stringr_1.4.0      dplyr_1.0.7
## [25] purrr_0.3.4        readr_2.1.1        tidyr_1.1.4        tibble_3.1.6
## [29] ggplot2_3.3.5      tidyverse_1.3.1
## loaded via a namespace (and not attached):
##  [1] VGAM_1.1-5           colorspace_2.0-2     ellipsis_0.3.2
##  [4] class_7.3-19         visdat_0.5.3         fs_1.5.2
##  [7] rstudioapi_0.13      listenv_0.8.0        MatrixModels_0.5-0
## [10] prodlim_2019.11.13   fansi_0.5.0          xml2_1.3.3
## [13] codetools_0.2-18     splines_4.1.2        Formula_1.2-4
## [16] jsonlite_1.7.2       pROC_1.18.0          mcmc_0.9-7
## [19] caret_6.0-90         rJava_1.0-6          broom_0.7.10
## [22] dbplyr_2.1.1         geepack_1.3-2        compiler_4.1.2
## [25] httr_1.4.2           backports_1.4.0      assertthat_0.2.1
## [28] survey_4.1-1         cli_3.1.0            quantreg_5.86
## [31] tools_4.1.2          coda_0.19-4          gtable_0.3.0
## [34] glue_1.5.1           reshape2_1.4.4       Rcpp_1.0.7
## [37] carData_3.0-4        cellranger_1.1.0     vctrs_0.3.8
## [40] Amelia_1.8.0         nlme_3.1-153         conquer_1.2.1
## [43] iterators_1.0.13     timeDate_3043.102    gower_0.2.2
## [46] globals_0.14.0       xlsxjars_0.6.1       rvest_1.0.2
## [49] lifecycle_1.0.1      future_1.23.0        ipred_0.9-12
## [52] miscTools_0.6-26     hms_1.1.1            parallel_4.1.2
## [55] SparseM_1.81         gridExtra_2.3        rpart_4.1-15
## [58] stringi_1.7.6        foreach_1.5.1        AER_1.2-9
## [61] boot_1.3-28          lava_1.6.10          matrixStats_0.61.0
## [64] rlang_0.4.12         pkgconfig_2.0.3      lattice_0.20-45
## [67] recipes_0.1.17       tidyselect_1.1.1     parallelly_1.29.0
## [70] plyr_1.8.6           magrittr_2.0.1       R6_2.5.1
## [73] generics_0.1.1       DBI_1.1.1            pillar_1.6.4
## [76] haven_2.4.3          foreign_0.8-81       withr_2.4.3
## [79] nnet_7.3-16          abind_1.4-5          future.apply_1.8.1
## [82] modelr_0.1.8         crayon_1.4.2         car_3.0-12
## [85] utf8_1.2.2           tzdb_0.2.0           maxLik_1.5-2
## [88] grid_4.1.2           data.table_1.14.2    ModelMetrics_1.2.2.2
## [91] reprex_2.0.1         digest_0.6.29        MCMCpack_1.6-0
## [94] stats4_4.1.2         MatchIt_4.3.2        munsell_0.5.0
## [97] mitools_2.4
# ----------------------------------------------------------------------
# load ST clean data
# ----------------------------------------------------------------------
## Uncomment and set working directory to replication folder
## setwd("folder")
rm(list = ls())
setwd("C:/Users/User/Desktop/EUI/2nd term/Replicating Research/10/BSS replication files/Weiss et al original files/doc_replication")
st_final <- readRDS("data/st_final.RDS")
# replace NA with 0
final_data <- st_final %>%
mutate(.,
daily_join = ifelse(is.na(daily_join), 0, daily_join),
join_binary = ifelse(daily_join >0,1,0),
week = lubridate::week(ymd(st_final$date)),
year = lubridate::year(ymd(st_final$date)),
month = lubridate::month(ymd(st_final$date))) %>%
mutate(year_fac = factor(year),
month_fac = factor(month),
week_fac = factor(week),
lcode_fac = factor(lcode))
# ----------------------------------------------------------------------
# Assign post deal date
# ----------------------------------------------------------------------
final_data <- final_data %>%
mutate(.,
post = ifelse(date > "2020-01-27",1,0))
# ----------------------------------------------------------------------
# Mixed and Arab only data
# ----------------------------------------------------------------------
# Create subset of non-jewish localities
final_data_nj <- final_data %>%
filter(.,
non_jewish == 1)
sum(final_data_nj$daily_join, na.rm=T)
nrow(final_data_nj)
# ----------------------------------------------------------------------
# Create descriptive stats table
# ----------------------------------------------------------------------
disc_table_nj <- dplyr::select(final_data_nj,
c(daily_join, join_binary,
triangle, ext_tri, pop_2018, age0_19sef_pcnt_08,
age65sef_pcnt_08, age85sef_pcnt_08,
AcadmCert_pcnt_08, Worked_08, HousingDens_avg_08,
Vehicle1up_pcnt_08, ChldBorn_avg_08))
stargazer(as.data.frame(disc_table_nj),
covariate.labels = c("Daily Join (Count)", "Daily Join (Binary)",
"Triangle", "Extended Triangle",  "Population 2018",
"Perc. Age 0-19", "Perc. Age 65+", "Perc. Age 85+",
"Perc. Academic", "Perc. Employed", "Housing Density",
"HH with Vehicle", "Average Children per Woman"),
title = "Descriptive Statistics - Non Jewish Localities",
label = "tab:dsc_st",
style = "qje",
omit.summary.stat = c("p25", "p75"),
notes = c("All variables following 'Perc. Age 0-19' are from the 2008 census."),
notes.append = T,
notes.align = "l",
out="tables/table_A9.tex")
# ----------------------------------------------------------------------
# Basic DiD
# ----------------------------------------------------------------------
# Event study
df_es <- final_data %>%
mutate(
event_time_day = as.numeric(date - as.Date("2020-01-28"))
) %>%
filter(event_time_day >= -60 & event_time_day <= 60)  # Keep only ±60 days
event_model_day <- feols(
join_binary ~ i(event_time_day, triangle, ref = 0) | lcode + event_time_day,
cluster = "lcode",
data = df_es
)
library(fixest)
## Event study
# days
df_es <- final_data %>%
mutate(
event_time_day = as.numeric(date - as.Date("2020-01-28"))
) %>%
filter(event_time_day >= -60 & event_time_day <= 60)  # Keep only ±60 days
event_model_day <- feols(
join_binary ~ i(event_time_day, triangle, ref = 0) | lcode + event_time_day,
cluster = "lcode",
data = df_es
)
iplot(event_model_day,
xlab = "Days since treatment",
ylab = "Effect on Joining Movement",
main = "Event Study: ±60 Days Around Treatment",
col = "firebrick3")
event_model_day <- feols(
join_binary ~ i(event_time_day, triangle, ref = -1) | lcode + event_time_day,
cluster = "lcode",
data = df_es
)
iplot(event_model_day,
xlab = "Days since treatment",
ylab = "Effect on Joining Movement",
main = "Event Study: ±60 Days Around Treatment",
col = "firebrick3")
png("plots/BSS_figure_3.png", width = 6.5, height = 4, units = "in", res = 300)
iplot(event_model_day,
xlab = "Days since treatment",
ylab = "Effect on Joining Movement",
main = "Event Study: ±60 Days Around Treatment",
col = "firebrick3")
dev.off()
# months
did_3 <- feols(
join_binary ~ i(month_seq, triangle, ref = -1) | lcode + month_seq,
cluster = "lcode",
data = df
)
# months
df <- final_data_nj %>%
arrange(year, month) %>%
mutate(
month_seq = (year - 2020) * 12 + (month - 1)
)
did_3 <- feols(
join_binary ~ i(month_seq, triangle, ref = -1) | lcode + month_seq,
cluster = "lcode",
data = df
)
iplot(did_3,
xlab = "Months since treatment (0 = Jan 2020)",
ylab = "Effect on Joining Movement",
main = "Event Study: Monthly Time Scale (Ref = Dec 2019)",
col = "firebrick3")
png("plots/BSS_figure_4.png", width = 6.5, height = 4, units = "in", res = 300)
iplot(did_3,
xlab = "Months since treatment (Jan 2020)",
ylab = "Effect on Joining Movement",
main = "Event Study: Monthly Time Scale Centered on Treatment",
col = "firebrick3")
dev.off()
library(lubridate)
library(scales)
final_data_nj$date <- as.Date(final_data_nj$date)
# after mid-2019, by week
pt <- final_data_nj %>%
filter(date > as.Date("2019-07-01")) %>%
mutate(week = floor_date(date, "week")) %>%
group_by(triangle, week) %>%
summarize(daily_join = sum(daily_join, na.rm = TRUE))
library(dplyr)
final_data_nj$date <- as.Date(final_data_nj$date)
# after mid-2019, by week
pt <- final_data_nj %>%
filter(date > as.Date("2019-07-01")) %>%
mutate(week = floor_date(date, "week")) %>%
group_by(triangle, week) %>%
summarize(daily_join = sum(daily_join, na.rm = TRUE))
# after mid-2019, by week
pt <- final_data_nj %>%
filter(date > as.Date("2019-07-01")) %>%
mutate(week = floor_date(date, "week")) %>%
group_by(triangle, week) %>%
dplyr::summarize(daily_join = sum(daily_join, na.rm = TRUE))
pt <- mutate(pt,
triangle = ifelse(triangle == 1, "Triangle", "Non-Triangle"))
pt$triangle <- factor(pt$triangle, levels = c("Triangle", "Non-Triangle"))
ggplot(pt, aes(x = week, y = daily_join, color = as.factor(triangle),
shape = as.factor(triangle))) +
geom_line(aes(color = as.factor(triangle))) +
geom_point(aes(color = as.factor(triangle))) +
theme(axis.text.x = element_text(size = 10)) +
geom_vline(xintercept = as.Date("2020-01-28"), linetype = "dotdash",
color = "grey", size = 0.6) +
annotate("text", x = as.Date("2019-09-01"), y = 1400, label = "PRE",
family = "Times", size = 3) +
annotate("text", x = as.Date("2020-04-01"), y = 1400, label = "POST",
family = "Times", size = 3) +
scale_color_manual(values = c("firebrick3", "dodgerblue3")) +
scale_shape_manual(values = c("triangle", "circle")) +
scale_y_continuous(labels = label_comma()) +
labs(color = "", shape = "",
x = "",
y = "Weekly Mobilization") +
theme(text = element_text(size = 12, family = "Times"),
legend.key = element_blank(),
panel.grid.major = element_blank(),
axis.text.x = element_text(size = 12),
plot.caption = element_text(size = 12, family = "Times", hjust = -.02),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"))
ggsave("plots/BSS_figure_5.png", width = 6.5, height = 4)
final_data_nj$join_binary_3plus <- ifelse(final_data_nj$daily_join >= 3, 1, 0)
final_data$join_binary_3plus <- ifelse(final_data$daily_join >= 3, 1, 0)
m1_thresh3 <- felm(join_binary_3plus ~ triangle * post + pop_2018 | 0 | 0 | lcode,
data = final_data_nj,
exactDOF = TRUE,
keepCX = TRUE)
summary(m1_thresh3)
## Table 2 model 2: main specification with time and locality fixed effects
m2 <- felm(join_binary_3plus ~ triangle * post | lcode + year + month + week | 0 | lcode,
data = final_data_nj,
exactDOF = TRUE,
keepCX = TRUE)
summary(m2)
## Table 2 model 3: alternative treatment with time and locality fixed effects
m3 <- final_data_nj %>%
mutate(.,
triangle = ext_tri) %>%
felm(join_binary_3plus ~ triangle * post | lcode + year + month + week | 0 | lcode,
data = .,
exactDOF = TRUE,
keepCX = TRUE)
summary(m3)
## Table 2 model 4: full sample with time and locality fixed effects
m4 <- felm(join_binary_3plus ~ triangle * post | lcode + year + month + week | 0 | lcode,
data = final_data,
exactDOF = TRUE,
keepCX = TRUE)
summary(m4)
# OUTPUT
se_list <- list(
summary(m1_thresh3)$coef[, "Cluster s.e."],
summary(m2)$coef[, "Cluster s.e."],
summary(m3)$coef[, "Cluster s.e."],
summary(m4)$coef[, "Cluster s.e."]
)
se_list <- lapply(se_list, function(x) { x[x < 1e-6] <- NA; x })
stargazer(m1_thresh3, m2, m3, m4,
out = "tables/table7_BSS.tex",
style = "qje",
title = "Deal of the Century Effect on Mobilization (3+ Subscriptions Threshold)",
keep = c("triangle", "post", "triangle:post"),
dep.var.labels = c("Mobilization (3+ joins)"),
covariate.labels = c("Triangle", "Post", "Triangle * Post"),
omit.stat = c("rsq", "f", "ser", "adj.rsq"),
ci = FALSE,
se = se_list,
omit.table.layout = "n",
column.sep.width = "-6pt",
star.cutoffs = NULL,
star.char = c("", "", ""),
label = "tab:main_did_st_thresh3",
add.lines = list(
c("Population Control", "Yes", "No", "No", "No"),
c("Week FE", "No", "Yes", "Yes", "Yes"),
c("Month FE", "No", "Yes", "Yes", "Yes"),
c("Year FE", "No", "Yes", "Yes", "Yes"),
c("Locality FE", "No", "Yes", "Yes", "Yes"),
c("Cluster", "Locality", "Locality", "Locality", "Locality"),
c("Sample", "Non-Jewish", "Non-Jewish", "Non-Jewish", "Full"),
c("Treatment", "10 Localities", "10 Localities", "16 Localities", "10 Localities"),
c("Pre-Register", "No", "No", "No", "No")
)
)
